function [alphas,betas,R2s, mats_yrs] = CSregs_ver2_nom(gesol, Phi_NxYz, stationary_prob_NxYz, Y_grid, FLAG_RUN_PARALLEL)
% This code implements the following version of Campbell-Shiller:
% hpxr(t->t+12,n) regressed on slope(t,y) = y(t,n) - y(t,1)

    if nargin<5; FLAG_RUN_PARALLEL = false; end
    
    mats_yrs = 2:5;
    
    alphas = zeros(size(mats_yrs));
    betas = zeros(size(mats_yrs));
    R2s = zeros(size(mats_yrs));

    Phi_NxYz_12mo = mpower2(Phi_NxYz,12);
    
    if FLAG_RUN_PARALLEL
        
        lhs0_list = cell(1,4);
        lhst_list = cell(1,4);
        rhs0_list = cell(1,4);
        
        inflation_process = gesol.inflation_process;
        
        for n = 2:5
            
            idx_buy = find(gesol.nom_bonds.maturities==12*n);
            idx_sell = find(gesol.nom_bonds.maturities==12*(n-1));
            idx_1y = find(gesol.nom_bonds.maturities==12);

            % -p(n)+p(1)
            lhs0.fun_Nxz = -log(gesol.nom_bonds.P_A(:,:,:,idx_buy)) + log(gesol.nom_bonds.P_A(:,:,:,idx_1y));
            lhs0.slope_Y = gesol.nom_bonds.P_C_a(idx_buy) - gesol.nom_bonds.P_C_a(idx_1y);
            lhs0.const = gesol.nom_bonds.P_B_a(idx_buy) - gesol.nom_bonds.P_B_a(idx_sell);
            lhs0.slope_X = gesol.nom_bonds.P_B_b(idx_buy) - gesol.nom_bonds.P_B_b(idx_sell);
            lhs0.slope_eps = gesol.nom_bonds.P_B_c(idx_buy) - gesol.nom_bonds.P_B_c(idx_sell);
            
            lhs0_list{n-1} = lhs0;

            % p(n-1)
            lhst.fun_Nxz = log(gesol.nom_bonds.P_A(:,:,:,idx_sell));
            lhst.slope_Y = -gesol.nom_bonds.P_C_a(idx_sell);
            lhst.const = -gesol.nom_bonds.P_B_a(idx_sell);
            lhst.slope_X = -gesol.nom_bonds.P_B_b(idx_sell);
            lhst.slope_eps = -gesol.nom_bonds.P_B_c(idx_sell);
            
            lhst_list{n-1} = lhst;

            % rhs = y(n)-y(1)
            rhs0.fun_Nxz = (-log(gesol.nom_bonds.P_A(:,:,:,idx_buy))/n ...
                + log(gesol.nom_bonds.P_A(:,:,:,idx_1y)));
            rhs0.slope_Y = (gesol.nom_bonds.P_C_a(idx_buy)/n - gesol.nom_bonds.P_C_a(idx_1y));
            rhs0.const = (gesol.nom_bonds.P_B_a(idx_buy)/n - gesol.nom_bonds.P_B_a(idx_1y));
            rhs0.slope_X = (gesol.nom_bonds.P_B_b(idx_buy)/n - gesol.nom_bonds.P_B_b(idx_1y));
            rhs0.slope_eps = (gesol.nom_bonds.P_B_c(idx_buy)/n - gesol.nom_bonds.P_B_c(idx_1y));
            
            rhs0_list{n-1} = rhs0;
            
        end
        
        parfor n = 2:5

            [alphas(n-1),betas(n-1),R2s(n-1)] = generic_nominal_reg(inflation_process, lhs0_list{n-1}, 12, lhst_list{n-1}, rhs0_list{n-1}, stationary_prob_NxYz, Phi_NxYz_12mo, Y_grid);

        end
        
    else
    
        for n = 2:5

            idx_buy = find(gesol.nom_bonds.maturities==12*n);
            idx_sell = find(gesol.nom_bonds.maturities==12*(n-1));
            idx_1y = find(gesol.nom_bonds.maturities==12);

            % -p(n)+p(1)
            lhs0.fun_Nxz = -log(gesol.nom_bonds.P_A(:,:,:,idx_buy)) + log(gesol.nom_bonds.P_A(:,:,:,idx_1y));
            lhs0.slope_Y = gesol.nom_bonds.P_C_a(idx_buy) - gesol.nom_bonds.P_C_a(idx_1y);
            lhs0.const = gesol.nom_bonds.P_B_a(idx_buy) - gesol.nom_bonds.P_B_a(idx_sell);
            lhs0.slope_X = gesol.nom_bonds.P_B_b(idx_buy) - gesol.nom_bonds.P_B_b(idx_sell);
            lhs0.slope_eps = gesol.nom_bonds.P_B_c(idx_buy) - gesol.nom_bonds.P_B_c(idx_sell);

            % p(n-1)
            lhst.fun_Nxz = log(gesol.nom_bonds.P_A(:,:,:,idx_sell));
            lhst.slope_Y = -gesol.nom_bonds.P_C_a(idx_sell);
            lhst.const = -gesol.nom_bonds.P_B_a(idx_sell);
            lhst.slope_X = -gesol.nom_bonds.P_B_b(idx_sell);
            lhst.slope_eps = -gesol.nom_bonds.P_B_c(idx_sell);

            % rhs = y(n)-y(1)
            rhs0.fun_Nxz = (-log(gesol.nom_bonds.P_A(:,:,:,idx_buy))/n ...
                + log(gesol.nom_bonds.P_A(:,:,:,idx_1y)));
            rhs0.slope_Y = (gesol.nom_bonds.P_C_a(idx_buy)/n - gesol.nom_bonds.P_C_a(idx_1y));
            rhs0.const = (gesol.nom_bonds.P_B_a(idx_buy)/n - gesol.nom_bonds.P_B_a(idx_1y));
            rhs0.slope_X = (gesol.nom_bonds.P_B_b(idx_buy)/n - gesol.nom_bonds.P_B_b(idx_1y));
            rhs0.slope_eps = (gesol.nom_bonds.P_B_c(idx_buy)/n - gesol.nom_bonds.P_B_c(idx_1y));

            [alphas(n-1),betas(n-1),R2s(n-1)] = generic_nominal_reg(gesol.inflation_process, lhs0, 12, lhst, rhs0, stationary_prob_NxYz, Phi_NxYz_12mo, Y_grid);

        end
    
    end

end

%% old save
% 
%     mats_yrs = 2:5;
%     
%     alphas = zeros(size(mats_yrs));
%     betas = zeros(size(mats_yrs));
%     R2s = zeros(size(mats_yrs));
% 
%     Phi_NxYz_12mo = mpower2(Phi_NxYz,12);
% 
%     [NN, Nx, NY, Nz] = size(stationary_prob_NxYz);
%     Y_NxYz = permute(repmat(Y_grid(:),[1 NN Nx Nz]),[2 3 1 4]);
% 
%     inflation_process = gesol.inflation_process;
%     mean_X = inflation_process.mu_pi/(1 - inflation_process.rho_pi);
%     var_X = (1 + 2*inflation_process.rho_pi*inflation_process.theta_pi + inflation_process.theta_pi^2)...
%         *(inflation_process.sigma_pi^2)/(1 - inflation_process.rho_pi^2);
%     var_e = inflation_process.sigma_pi^2;
%     cov_Xe = var_e;
%     
%     for i = 1:numel(mats_yrs)
%         
%         mat_yr = mats_yrs(i);
%         mat_mth = mat_yr*12;
%         
%         idx_1y = find(gesol.nom_bonds.maturities==12);
%         idx_ny = find(gesol.nom_bonds.maturities==12*mat_yr);
%         idx_nless1y = find(gesol.nom_bonds.maturities==12*(mat_yr-1));
%         
%         if isempty(idx_1y) || isempty(idx_ny) || isempty(idx_nless1y)
%             error('Some yields are not available');
%         end
%         
%         % compute coefficients for RHS: slope = y(n) - y(1)
%         rhs_Nxz = -log(squeeze(gesol.nom_bonds.P_A(:,:,:,idx_ny)))/mat_yr ...
%             + log(squeeze(gesol.nom_bonds.P_A(:,:,:,idx_1y)));
% 
%         rhs_coeff_Y = gesol.nom_bonds.P_C_a(idx_ny)/mat_yr ...
%             - gesol.nom_bonds.P_C_a(idx_1y);
%         rhs_coeff_X = -gesol.nom_bonds.P_B_b(idx_ny)/mat_yr ...
%             - gesol.nom_bonds.P_B_b(idx_1y);
%         rhs_coeff_eps = -gesol.nom_bonds.P_B_c(idx_ny)/mat_yr ...
%             - gesol.nom_bonds.P_B_c(idx_1y);
% 
%         rhs_NxYz = permute(repmat(rhs_Nxz,[1 1 1 NY]),[1 2 4 3]) + rhs_coeff_Y*Y_NxYz;
%         
%         % LHS = p(n-1,t+1) - p(n,t) + p(1,t)
%         % LHS0 = - p(n,t) + p(1,t), LHS12 at h=12 is p(n-1,t+1)
%         
%         lhs0_Nxz = -log(squeeze(gesol.nom_bonds.P_A(:,:,:,idx_ny))) ...
%             + log(squeeze(gesol.nom_bonds.P_A(:,:,:,idx_1y)));
%         lhs0_coeff_Y = gesol.nom_bonds.P_C_a(idx_ny) - gesol.nom_bonds.P_C_a(idx_1y);
%         lhs0_coeff_X = gesol.nom_bonds.P_B_b(idx_ny) - gesol.nom_bonds.P_B_b(idx_1y);
%         lhs0_coeff_eps = gesol.nom_bonds.P_B_c(idx_ny) - gesol.nom_bonds.P_B_c(idx_1y);
% 
%         lhs0_NxYz = permute(repmat(lhs0_Nxz,[1 1 1 NY]),[1 2 4 3]) + lhs0_coeff_Y*Y_NxYz;
% 
%         lhs12_Nxz = log(squeeze(gesol.nom_bonds.P_A(:,:,:,idx_nless1y)));
%         lhs12_coeff_Y = -gesol.nom_bonds.P_C_a(idx_nless1y);
%         lhs12_coeff_X = -gesol.nom_bonds.P_B_b(idx_nless1y);
%         lhs12_coeff_eps = -gesol.nom_bonds.P_B_c(idx_nless1y);
% 
%         lhs12_NxYz = permute(repmat(lhs12_Nxz,[1 1 1 NY]),[1 2 4 3]) + lhs12_coeff_Y*Y_NxYz;
%         
%         % reg coeffs
%         mean_rhs_NxYz = sum(stationary_prob_NxYz(:).*rhs_NxYz(:));
%         var_rhs_NxYz = sum(stationary_prob_NxYz(:).*(rhs_NxYz(:).^2)) - mean_rhs_NxYz^2;
%         mean_rhs_Xe = rhs_coeff_X*mean_X;
%         var_rhs_Xe = (rhs_coeff_X^2)*var_X + 2*rhs_coeff_X*rhs_coeff_eps*cov_Xe ...
%             + (rhs_coeff_eps^2)*var_e;
%         mean_rhs = mean_rhs_NxYz + mean_rhs_Xe;
%         var_rhs = var_rhs_NxYz + var_rhs_Xe;
% 
%         mean_lhs0_NxYz = sum(stationary_prob_NxYz(:).*lhs0_NxYz(:));
%         mean_lhs12_NxYz = sum(stationary_prob_NxYz(:).*lhs12_NxYz(:));
%         mean_lhs0_Xe = lhs0_coeff_X*mean_X;
%         mean_lhs12_Xe = lhs12_coeff_X*mean_X;
%         mean_lhs = mean_lhs0_NxYz + mean_lhs12_NxYz + mean_lhs0_Xe + mean_lhs12_Xe;
% 
%         var_lhs0_NxYz = sum(stationary_prob_NxYz(:).*(lhs0_NxYz(:).^2)) - mean_lhs0_NxYz^2;
%         var_lhs12_NxYz = sum(stationary_prob_NxYz(:).*(lhs12_NxYz(:).^2)) - mean_lhs12_NxYz^2;
%         cov_lhs0_lhs12_NxYz = sum(stationary_prob_NxYz(:).*lhs0_NxYz(:).*(Phi_NxYz_12mo*lhs12_NxYz(:))) ...
%             - mean_lhs0_NxYz*mean_lhs12_NxYz;
%         var_lhs_NxYz = var_lhs0_NxYz + var_lhs12_NxYz + 2*cov_lhs0_lhs12_NxYz;
% 
%         var_lhs0_Xe = (lhs0_coeff_X^2)*var_X + (lhs0_coeff_eps^2)*var_e ...
%             + 2*lhs0_coeff_X*lhs0_coeff_eps*cov_Xe;
%         var_lhs12_Xe = (lhs12_coeff_X^2)*var_X + (lhs12_coeff_eps^2)*var_e ...
%             + 2*lhs12_coeff_X*lhs12_coeff_eps*cov_Xe;
%         cov_lhs0_lhs12_Xe = lhs12_coeff_X*lhs0_coeff_X*((inflation_process.rho_pi^12)*var_X ...
%             + (inflation_process.rho_pi^11)*inflation_process.theta_pi*cov_Xe) ...
%             + lhs12_coeff_X*lhs0_coeff_eps*((inflation_process.rho_pi^12)*cov_Xe ...
%             + (inflation_process.rho_pi^11)*inflation_process.theta_pi*var_e);
%         var_lhs_Xe = var_lhs0_Xe + var_lhs12_Xe + 2*cov_lhs0_lhs12_Xe;
%         var_lhs = var_lhs_NxYz + var_lhs_Xe;
% 
%         cov_rhs_lhs0_Xe = lhs0_coeff_X*rhs_coeff_X*var_X + (lhs0_coeff_X*rhs_coeff_eps ...
%             + lhs0_coeff_eps*rhs_coeff_X)*cov_Xe + lhs0_coeff_eps*rhs_coeff_eps*var_e;
%         cov_rhs_lhs12_Xe = lhs12_coeff_X*rhs_coeff_X*((inflation_process.rho_pi^12)*var_X ...
%             + (inflation_process.rho_pi^11)*inflation_process.theta_pi*cov_Xe) ...
%             + lhs12_coeff_X*rhs_coeff_eps*((inflation_process.rho_pi^12)*cov_Xe ...
%             + (inflation_process.rho_pi^11)*inflation_process.theta_pi*var_e);
% 
%         cov_rhs_lhs0_NxYz = sum(stationary_prob_NxYz(:).*rhs_NxYz(:).*lhs0_NxYz(:)) ...
%             - mean_rhs_NxYz*mean_lhs0_NxYz;
%         cov_rhs_lhs12_NxYz = sum(stationary_prob_NxYz(:).*rhs_NxYz(:).*(Phi_NxYz_12mo*lhs12_NxYz(:))) ...
%             - mean_rhs_NxYz*mean_lhs12_NxYz;
% 
%         cov_lhs_rhs = cov_rhs_lhs0_NxYz + cov_rhs_lhs12_NxYz ...
%             + cov_rhs_lhs0_Xe + cov_rhs_lhs12_Xe;
%         
%         betas(i) = cov_lhs_rhs/var_rhs;
%         alphas(i) = mean_lhs - betas(i)*mean_rhs;
%         R2s(i) = (betas(i)^2)*var_rhs/var_lhs;
%         
%     end